WT RPM Organoids Allografts Notebook

Modifying Original file (Fig2_ExtFig4_WT_RPM_Organoids_Allografts_Final.R) into an executable RMarkdown Notebook

Related to:

• Fig 1i-n
• Extended Data Fig 4a-f

suppressPackageStartupMessages({
    library(viridis)
    library(ggplot2)
    library(scater)
    library(Seurat)
    library(ggpubr)
})

Define colors for figures

# Define color scheme and plot UMAP as in Fig. 2e #
my_colors <- c(
  "#E41A1C", # strong red
  "#377EB8", # medium blue
  "#4DAF4A", # green
  "#984EA3", # purple
  "#FF7F00", # orange
  "#FFFF33", # yellow
  "#A65628", # brown
  "#e7298a", # pink
  "#666666", # grey
  "lavender", # teal
  "#FC8D62", # salmon
  "#8DA0CB", # soft blue
  "#E78AC3", # soft pink (different from 8)
  "#A6D854", # light green (but yellowish tint, not green)
  "#FFD92F", # lemon yellow
  "#E5C494", # light brown
  "#B3B3B3", # light grey
  "#1B9E77", # deep teal
  "#D95F02", # dark orange
  "#7570B3", # strong purple
  "turquoise"  # olive green (NOT same green as before)
)

colors_figED4a <- c('deepskyblue4', '#FC8D62','#2ca02c', 'brown1', 'purple', '#A0522D',
                    'pink2',  '#a1c299',  '#0aa6d8', 'lightblue','salmon1','lightgreen')

sample_cols<-c("orange","#D8BFD8","darkorchid4", # deep teal
               "turquoise", "#1B9E77" # dark orange
)

pheno_col<-c("brown2","darkorchid4","dodgerblue","#66A61E","orange","turquoise4")

Load RPM Organoids data

RPM_Orgs <- readRDS("../data/05_2025_RPM_Orgs_Only_ExtFig4a-c.rds")

Load CellTagged Organoid Allografts data

RPM_Orgs_Allo <- readRDS("../data/05_2025_RPM_Orgs_Allo_Fig2d.rds")

Load RPM allograft-only data

RPM_Allo <- readRDS("../data/05_2025_RPM_AllograftOnly_Fig2e.rds")

Load RPM allografts as SingleCellExperiment

sce <- readRDS('../data/RPM_Allo_only_sce.rds')

UMAP Plots

Fig 1i

DimPlot(RPM_Orgs_Allo,  group.by='Sample',cols=sample_cols, reduction='umap',
        label=FALSE, label.size=6, shuffle=TRUE, pt.size=0.01) & NoAxes()

FeaturePlot(RPM_Orgs_Allo, features = c("NE_Consensus1"), 
            pt.size=0.2, reduction='umap', order=TRUE) + 
    ggplot2::scale_color_gradientn(colors=rocket(10, direction=-1)) & NoAxes()

FeaturePlot(RPM_Orgs_Allo, features = c("Basal_Consensus1"), 
            pt.size=0.2, reduction='umap',order=TRUE) + 
    ggplot2::scale_color_gradientn(colors=rocket(10, direction=-1)) & NoAxes()

Fig 1j

Violin plots of expression of TFs by Leiden cluster

genes <- c("Ascl1","Neurod1","Pou2f3","Atoh1","Yap1","Trp63")

# Create violin plots with Kruskal-Wallis test
plots <- lapply(genes, function(gene) {
  p <- VlnPlot(RPM_Orgs_Allo,
               features = gene,
               group.by = "leiden_scVI_1.3",
               cols = my_colors,
               alpha = 0.5) +  # smaller dots
    ggpubr::stat_compare_means(method = "kruskal.test", 
                       label = "p.format", 
                       label.x = 2,size = 5) +
    ggtitle("") +
    theme(plot.title = element_text(face = "italic"), legend.position = "none",  # Remove legend
          axis.title.x = element_blank(),axis.title.y = element_text(size = 20),
          axis.text.y = element_text(size = 20), axis.text.x = element_text(size = 14)) +  # Remove x-axis label
    labs(y = "Expression")  # Change y-axis label to "Expression"
  return(p)
})

# Arrange plots
patchwork::wrap_plots(plots, ncol = 6)

UMAP Plots by Sample and Cluster

Extended Data Fig 4a

DimPlot(RPM_Orgs, group.by='Genotype', cols=c("darkorchid4","orange"), reduction='umap',
        label=FALSE, label.size=6, shuffle=TRUE) & NoAxes()

DimPlot(RPM_Orgs, group.by='cluster', cols=colors_figED4a, reduction='umap',
        label=TRUE, label.size=10) & NoAxes()

Extended Data Fig 4c

DimPlot(RPM_Orgs, reduction = "umap", group.by = "Phase",
        shuffle=TRUE, label=FALSE, pt.size=.05, 
        cols=c("indianred3","green3","royalblue4")) + NoAxes()

Ext. Data Fig. 4d

DimPlot(RPM_Orgs_Allo, reduction = "umap", group.by = "Phase",
        shuffle=TRUE, label=FALSE, pt.size=.05, 
        cols=c("indianred3","green3","royalblue4")) + NoAxes()

# Plot barplot for Ext. Data Fig. 4d #
##### What % of cells per phase in each group? ####
x <- table(Idents(RPM_Orgs_Allo),RPM_Orgs_Allo@meta.data$Phase)
proportions <- as.data.frame(100*prop.table(x, margin = 1))

colnames(proportions)<-c("Cluster", "Sample", "Frequency")

# ggbarplot(proportions, x="Sample", y="Frequency", fill = "Sample", group = "Sample", ylab = "Frequency (percent)", xlab="Phase", palette =c("indianred3","green3","royalblue4"))+ theme_bw()+ facet_wrap(facets = "Cluster", scales="free_y", ncol =4)+ theme(axis.text.y = element_text(size=12)) +rotate_x_text(angle = 45)

# Stacked
p<-ggplot(proportions, aes(fill=Sample, y=Frequency, x=Cluster)) +
  geom_bar(position="stack", stat="identity")

p + scale_fill_manual(values=c("indianred3","green3","royalblue4")) + 
    theme_bw() + 
    theme(axis.text.y = element_text(size=20), axis.text.x=element_text(size=20,angle=45, hjust = 1), 
          axis.title.x = element_text(size=20), axis.title.y = element_text(size=20), 
          legend.text = element_text(size=20), legend.title = element_text(size=20)) +
    labs(x=NULL, y="Frequency (%)", fill=NULL)

Fig 1j

DimPlot(RPM_Allo, group.by='leiden_scVI_1.3', cols=my_colors,
        reduction='umap', label=TRUE, label.size=7) & NoAxes() + theme(legend.position="none")

UMAP by Fate

Fig. 1k

pheno_col<-c("brown2","darkorchid4","dodgerblue","#66A61E","orange","turquoise4")

DimPlot(RPM_Allo, group.by=c("Pheno"), cols=pheno_col, shuffle=TRUE, pt.size=0.6)+NoAxes()

Fig 1l

# Fig 2g plots #
FeaturePlot(RPM_Allo, features = c("NE_spearman"), pt.size=0.2, 
            reduction='umap', order=TRUE) + 
    scale_color_viridis(option="rocket",direction=-1) & NoAxes()

Violin plots of Archetypes by leiden cluster

Ext. Data Fig. 4f

VlnPlot(
  RPM_Allo,
  features = c("A_Archetype1","A2_Archetype1","N_Archetype1","P_Archetype1"),
  group.by = "leiden_scVI_1.3",
  cols = my_colors,
  alpha = 0.7,
  ncol = 2
) & 
  theme(
    axis.text.x = element_text(size = 8),                # smaller x-axis labels
    plot.title = element_text(face = "bold", size=0),          # italicize titles
    axis.title.x = element_text(size = 24),              # optional size tweaks
    axis.title.y = element_text(size = 24)
  ) &
  labs(
    y = "Expression",
    x = "Cluster"
  )

Plot Signature Scores by Phenotype

phenotypes <- as.character(unique(RPM_Allo$Pheno))
pheno_pairs_list <- combn(phenotypes, 2, simplify = FALSE)

Define function to plot violin plots by phenotype

plotVlnByPhen <- function(feature, yl=c(-0.1,0.3), focus=NA) {
    if(!is.na(focus)) {
        pheno_pairs_list <- pheno_pairs_list[sapply(pheno_pairs_list, function(x) focus %in% x)]
    }
    
    scater::plotColData(sce, x = "Pheno", y = feature, colour_by = "Pheno") + 
        scale_discrete_manual(aesthetics = c("colour", "fill"), values=pheno_col) +
        theme(axis.title.y = element_text(size = 24), axis.text.y = element_text(size = 16),
              axis.text.x = element_blank(),   # rotate x-axis labels
              plot.title = element_blank(),    # remove plot title
              legend.position = "none"         # remove legend
        ) + labs(x = "",                     # custom x-axis title
                 y = "Signature score"        # custom y-axis title
        ) + ylim(yl[1], yl[2]) + 
        geom_boxplot(fill=pheno_col, alpha=1/5, position = position_dodge(width = .2),
                     size=0.2, color="black", notch=TRUE, notchwidth=0.3, outlier.shape = 2, outlier.colour=NA) #   + 
        # stat_summary(fun = mean,
        #              geom = "point",
        #              shape = 18,
        #              size = 2,
        #              color = "red",
        #              position = position_dodge(width = 0.5)) + 
        # ggpubr::stat_compare_means(method = "wilcox.test", label = "p.signif", size=4, 
        #                            vjust = 0.5,
        #                            hide.ns = TRUE, comparisons = pheno_pairs_list)
}

Fig 1l

NE correlation

plotVlnByPhen("NE_spearman", yl=c(-.75, .8), focus="NE")

Fig 1m

SCLC cell type consensus signature scores * NE cell * Tuft cell * Basal cell

plotVlnByPhen("NE_Consensus1", yl=c(-0.5, 1))

plotVlnByPhen("Tuft_Consensus1", yl=c(-0.1, 0.5))

plotVlnByPhen("Basal_Consensus1", yl=c(-0.1, 1.1))

Fig 1n

TF target gene scores * ASCL1 targets * NEUROD1 targets * ATOH1 targets * POU2F3 targets * YAP1 activity score

plotVlnByPhen("ASCL1_Targets1", yl=c(-0.1,0.6))

plotVlnByPhen("NEUROD1_Targets1", yl=c(-0.1,0.25))

plotVlnByPhen("ATOH1_Targets1", yl=c(-0.1,.3))

plotVlnByPhen("POU2F3_Targets1", yl=c(-0.1, 0.22))

plotVlnByPhen("YAP1_Activity_Score1", yl=c(-0.1, 1))